Mereghetti etal. BMC Biophysics 2014, 7:4 
http://www.biomedcentral.eom/2046-1682/7/4 



RESEARCH ARTICLE Open Access 



Long range Debye-Huckel correction for 
computation of grid-based electrostatic forces 
between biomacromolecules 

Paolo Mereghetti 1 ' 2 , Michael Martinez 1 and Rebecca C Wade 13 " 



Abstract 

Background: Brownian dynamics (BD) simulations can be used to study very large molecular systems, such as 
models of the intracellular environment, using atomic-detail structures. Such simulations require strategies to contain 
the computational costs, especially for the computation of interaction forces and energies. A common approach is to 
compute interaction forces between macromolecules by precomputing their interaction potentials on 
three-dimensional discretized grids. For long-range interactions, such as electrostatics, grid-based methods are 
subject to finite size errors. We describe here the implementation of a Debye-Huckel correction to the grid-based 
electrostatic potential used in the SDA BD simulation software that was applied to simulate solutions of bovine serum 
albumin and of hen egg white lysozyme. 

Results: We found that the inclusion of the long-range electrostatic correction increased the accuracy of both the 
protein-protein interaction profiles and the protein diffusion coefficients at low ionic strength. 

Conclusions: An advantage of this method is the low additional computational cost required to treat long-range 
electrostatic interactions in large biomacromolecular systems. Moreover, the implementation described here for BD 
simulations of protein solutions can also be applied in implicit solvent molecular dynamics simulations that make use 
of gridded interaction potentials. 

Keywords: Continuum solvent electrostatics, Ionic strength, Debye-Huckel , Poisson-Boltzmann equation, Brownian 
dynamics simulation, Protein diffusion, Discretization grid, Finite difference, Second virial coefficient, Small angle 
scattering intensity 



Background 

Simulations of concentrated solutions of macromolecules 
such as those designed to mimic the intracellular envi- 
ronment, are becoming feasible due to improvements 
in computational power and simulation methods [1-5]. 
Given that even for simulating a small volume of a protein 
solution, several hundreds of proteins have to be taken 
into account, coarse grained methods, which neglect 
atomic details, e.g. by treating each protein as a sphere, are 
often applied [6] . 
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However, to understand the effects of differences in 
protein sequence or point mutations from simulations 
requires a more detailed level of modeling. Explicit inclu- 
sion of atomic detail can be computationally demanding 
and therefore, approximations and calculation strategies 
are required to make the simulations feasible. A com- 
monly employed approach is to retain atomic detail for 
the macromolecules while treating them as rigid bodies 
in continuum solvent. Apart from restricting the number 
of degrees of freedom considered in the simulations, 
this treatment permits interaction forces between macro- 
molecules to be computed efficiently by precomputation 
of their interaction potentials on three-dimensional dis- 
cretized grids. Thus, during the simulations, forces can be 
computed by considering the interactions of each atom 
of each macromolecule with the interaction potential 
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grids of the other macromolecules. Grid formalisms 
for intermolecular interactions are extensively used for 
macromolecular docking methodologies [7,8], binding 
site determination [9], as well as in structure determi- 
nation from electron microscopy maps [10,11]. A major 
drawback of gridded potentials is, however, the occur- 
rence of finite size problems [3] . To minimize truncation 
errors in computing energies or forces, the interaction 
potential must be small at the edges of a grid. For molec- 
ular electrostatic potentials, the long-range nature of the 
Coulombic interaction, especially at low salt concentra- 
tion or for highly charged macromolecules, means that 
very large grids are often required. For example, at 5 mM 
ionic strength, the Debye length of the solution is 43 A. 
For a small globular protein with a radius of 20 A and a net 
charge of +10e, the electrostatic grid dimensions should 
be at least 200 x 200 x 200 A to obtain an electrostatic 
potential of ^0.1 kcal/mol/e at the grid edges. Assuming a 
grid spacing of 1 A, the grid must have at least 201 x 201 x 
201 points. This grid size is not a problem when a sin- 
gle small protein is considered but becomes an issue when 
simulating a periodic box containing several hundreds or 
thousands of proteins in solution. The grid size can also 
be a problem for memory usage in calculations for one or 
a few large macromolecules. 

One solution to this problem is to use multiple focused 
grids with different grid spacings centered on each macro- 
molecule: a detailed potential grid with a small grid 
spacing for representing the electrostatic potential at 
short-range and a coarse grid with a larger grid spac- 
ing for the long-range part [1]. Another solution, which 
will be described in this paper, is to exploit the fact 
that beyond a certain distance from the surface of the 
macromolecule, the electrostatic potential becomes cen- 
trosymmetric. Thus, a cubic gridded potential is used for 
the short-range part of the electrostatic potential up to 
a defined distance threshold and a continuous screened 
Coulomb potential is used beyond this distance. The dis- 
tance threshold corresponds to the radius of the largest 
sphere enclosed by the grid. 

We have recently developed a Brownian dynamics (BD) 
method for simulating many macromolecules (10 2 — 10 3 ) 
described as atomically detailed rigid bodies in a contin- 
uum solvent in a periodic box [3]. The model used is based 
on that originally developed for the simulation of the 
diffusional association of two proteins and implemented 
in the SDA (Simulation of Diffusional Association) soft- 
ware [8] . For the simulation of many proteins, this method 
gives results in good agreement with experimental trans- 
lational and rotational diffusion coefficients and small 
angle scattering structure factors for dilute [3] as well 
as concentrated protein solutions [12]. In this approach, 
intermolecular forces are computed as the sum of elec- 
trostatic interaction, electrostatic desolvation, non-polar 



desolvation and soft-core repulsion terms [3,8]. For com- 
putational efficiency, all these terms are precomputed on 
grids for each macromolecular solute before carrying out 
the BD simulations. To overcome errors due to the finite 
size of the electrostatic grids, we here describe the imple- 
mentation of a long-range electrostatic correction to the 
model for interaction forces used in our BD simulations. 
The purpose of this correction is to improve the accu- 
racy of the computed inter-protein forces and extend the 
applicability of the approach to highly charged proteins 
and low ionic strength conditions. For validation, we per- 
formed BD simulations of bovine serum albumin (BSA) 
and hen egg white lysozyme (HEWL) with and without 
the long-range electrostatic correction and compared the 
results to experimentally determined small angle scatter- 
ing structure factors and self-diffusion coefficients. The 
same methodology described here for the implementa- 
tion of the long-range Debye-Huckel correction, should 
also be applicable in implicit solvent molecular dynam- 
ics simulations that make use of gridded interaction 
potentials [13-16]. 

Methods 

Brownian dynamics (BD) is a simulation method that 
employs a mesoscopic model in which the solvent is 
treated as a continuum and the solutes are modelled as 
discrete entities at a level of detail appropriate for the 
problem being studied. BD thus takes advantage of the 
large separation in timescale between the fast solvent 
motion and the slower motion of solute particles (poly- 
mers or colloids) which make it possible to treat the 
solvent implicitly. Furthermore, internal solute degrees 
of freedom are often neglected and macromolecules 
are treated as rigid bodies interacting by direct inter- 
actions (electrostatic, van der Waals, non-polar) and 
solvent-mediated (hydrodynamic) interactions. Due to 
these simplifications, BD can be used to study larger 
biomacromolecular systems on longer time-scales than is 
possible with classical atomic-detail molecular dynamics 
simulations. 

Translational motion is propagated according to the 
following Equation [17]: 

r.'(fl) = r/fo) + T -——^—At + -£F ; -(ft>)Af + R; (1) 

*y 3tj(to) kT 

where ri is the position of the center of geometry of the 
solute i and At = (t\ — to) is the timestep. 

The effect of the solvent is described by a random 
displacement, R;, which mimics the collision of solute i 
with the solvent molecules and is defined by a Gaussian 
distribution with mean (R/) = 0 and covariance (R/R/) = 
2D; A£. From the latter, it follows that the stochastic 
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displacement is proportional to the square root of the 
translational diffusion tensor, DL The second term on 
the r.h.s. of Equation 1, the divergence of the diffusion 
tensor, describes the hydrodynamic drift of the solute 
towards regions of high mobility. The force acting on 
solute i results from the sum of the forces acting on 
solutes j at time to, Fy(£n)> coupled with the diffusion 
tensor. 

We employ a simplified treatment of hydrodynamic 
interactions to avoid the computationally expensive 
Cholesky factorization required to calculate the square 
root of the diffusion matrix. A mean field approach is 
used where D\- is replaced by a volume fraction-dependent 
diffusion coefficient, D tshort ((f)i), and Equation 1 simplifies 
to [12] 



Btshortffa) 

tiih) = n(to) + r^Fifo) At + R; 

kT 



(2) 



We define the local volume, V;, as the volume of the 
sphere of radius R cut centered on solute i. The local 
volume fraction 0; for the solute i is obtained by divid- 
ing the sum of the volumes of the solutes within R cut 
by the local volume V; [18]. The volume of a pro- 
tein, v, is computed by approximating the protein as 
a sphere having a radius equal to the hydrodynamic 
radius (cr stokes ) estimated using HYDROPRO [19]. The 
cutoff for the local volume, R cut , is set to four times 
the side of the largest interaction grid of the cen- 
tral solute. For a small simulation box, this cutoff was 
rescaled to a value equal to half of the simulation box 
size. A solute j is totally included in the local volume 
when the center-to-center distance between the cen- 
tral solute i and solute j is less than R cut — cr- tokes . 
When a solute k is only partially included within R cut , 



that is, when R cut - af kes < d ik < R cut + a k 
we account for that portion of solute volume derived 
by the sphere-sphere intersection. The volume frac- 
tion dependent short-time translational diffusion coeffi- 
cient (D tshort (00) is then obtained using the Tokuyama 
model [20-22], derived for a concentrated hard-sphere 
suspension of particles interacting with both direct and 
hydrodynamic interactions. An equation analogous to 
Equation 2 is used for the rotational motion [12], with 
the volume fraction dependent short-time rotational 
diffusion coefficient obtained using the model derived 
by Cichocki et al. which includes lubrication forces as 
well as two- and three-body expansions of the mobility 
functions [23]. 

The forces, Fi, are computed as finite-difference deriva- 
tives of the pairwise free energies of interaction between 
the solutes as described in the next section. 



stokes 



Interaction energies and forces 

For each pair of macromolecules, the interaction free 
energy, AG 1-2 , is defined as: 

AGl_2 = ^I>e/i(**2)-tt2 

h 

+ i ^ ®el 2 ( r /i) " [electrostatic interaction] 

h 

h 

+ ^ ®edesolv 2 ( r /i ) ' #/i [electrostatic desolvation] 

h 

+ y ] ^npdesolvi ( r W2 ) ' SASA m2 

+ ^ ®n P desolv 2 ( r «i ) • SASA ni [non-polar desolvation] 

m 

+ y ] Esoftcorei ( r W2 ) 
m 2 

+ X E soficore2 (r ni ) [soft-core repulsion] 

n\ 

(3) 

A detailed description and parameterization of 
Equation 3 can be found in Refs. [3,24]. Briefly, the first 
two terms in Equation 3 are the interaction energies 
of the charges of one macromolecule (qi 2 or qj x ) with 
the electrostatic potential of the other macromolecule 
or O e / 2 ). Charges were assigned using the effective 
charge approximation [25] . The third and fourth terms of 
Equation 3 represent the electrostatic desolvation energy 
arising from the introduction of the low dielectric cavity 
of one macromolecule in the presence of the charges of 
the other [25,26]. The desolvation energy is computed 
as the interaction of the charges of one macromolecule 
(qi 2 or qj x ) with the electrostatic desolvation potential 
of the other macromolecule {® e desolvi or ® e desolv 2 ) [ 26 1> 
with parameterization as in Ref. [24]. The fifth and sixth 
terms in Equation 3 correspond to the non-polar inter- 
actions due to the burial of the solvent accessible surface 
areas (SASAs) of the surface atoms The last two terms 
of Equation 3 describe the soft-core repulsive potential 
introduced to avoid overlaps. The soft-core potential is 
modelled using an inverse power function. The smooth- 
ness of the soft-core potential allows abrupt changes in 
the forces at close contact to be avoided. In Equation 3, 
r specifies the atomic coordinates. For computational 
efficiency, all interaction potentials, O, are mapped onto 
grids centered on each of the macromolecules. 

This formalism implies a truncation of the electro- 
static potential in the grid-charge formalism due to the 
finite extent of the grids. To alleviate this problem, we 
here introduce an analytical long-range correction to 
the electrostatic interaction term that makes use of the 
assumption that beyond the electrostatic grid boundaries, 
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a macromolecule can be treated as a Debye-Huckel 
sphere. 

According to the Debye-Huckel theory of dilute elec- 
trolyte solutions, all ions in the solvent are treated as 
point charges while each pair of solutes are treated 
as spheres with radii au aj and net charges z/e/, Zjei, 
where e\ is the elementary charge. Then, the poten- 
tial of mean force between a pair of solute molecules 
is 



w(r) = +oo 



2 



K{r-a) 



4?7T€Q€ r r{\+Ka) 



r < a 
r >a 



(4) 



Where r is the center-to-center distance and w(r) is the 
potential of mean force. For an isotropic potential, the 
corresponding equation is 



B22 



-Ul 



'< B t - l 



4?nr 2 dr 



(6) 



Small angle scattering intensity 

To assess the correctness of the interaction potentials, we 
compared experimental and computed small angle scat- 
tering intensities. Scattering intensities were computed 
from the simulations using [30] 



where 60 is the vacuum permittivity, e r is the relative per- 
mittivity of the solvent, a = at + aj, and k is the inverse of 
the Debye length, and is proportional to the ionic strength 



I(q) = yn p (Ap) 2 v 2 P(q)S(q) 



(« 2 = &^^> 



As shown in Equation 3, to compute the electrostatic 
interaction between a pair of macromolecules, the elec- 
trostatic potential of macromolecule 1 is multiplied by 
the effective charges of the second macromolecule. Due 
to the finite size of the grid, when the second macro- 
molecule is on the border of the electrostatic potential 
grid of macromolecule 1, only a fraction of the effec- 
tive charges on macromolecule 2 are taken into account 
for computing the electrostatic interaction. An isotropic 
distance cut-off from the center of macromolecule 1 
is used in computing this interaction, so that if the 
effective charge is beyond this distance cutoff, its elec- 
trostatic interaction is not computed. The spherical cut- 
off is assigned on the assumption that the electrostatic 
potential becomes centrosymmetric at the grid edges 
and therefore a switch to the analytical Debye-Huckel 
potential can be made beyond the cutoff. The appli- 
cation of the Debye-Huckel potential reduces the dis- 
continuity in the energy and forces at the grid cut-off 
distance. 

Second osmotic virial coefficients 

Osmotic virial coefficients are coefficients in the virial 
expansion of the state equation and they reflect deviations 
from ideal behaviour due to the presence of interactions. 
For simple cases, they can be obtained analytically. For 
this reason, they are commonly used to assess force field 
accuracy [1,3,27,28]. 

From classical statistical mechanics, the second osmotic 
virial coefficient can be obtained from [29] 



B22 



Mr) 
~k R T 



- 1 



d£2 



(5) 



(7) 



where 7 is a factor related to instrument effects, n p = 
N/V is the protein concentration expressed as number 
density (N is the number of particles and V the total 
volume of the solution), Ap is the electron density con- 
trast between the scattering particle and the solvent, and 
v is the particle volume. P(q) is the form factor normal- 
ized such that P(0) = 1, S(q) is the structure factor 
and q is the scattering vector. The pre-factor y(Ap) 2 v 2 
can be obtained in experiments and then the normalized 
scattering intensity is expressed as 



= P(q)S(q), whereA = y(Ap)V 



(8) 



We computed the form factor for BSA using the analyti- 
cal expression for the orientationally averaged form factor 
of an oblate ellipsoid with radii a and b where a is the 
semi-axis of revolution [31,32]. Following ref. [32], we set 
a = 17.5 A and b = 47.4 A. 

The structure factor, S(q), was computed by Fourier 
transformation of the radial distribution function, 
g(r) [33] as follows 



S(q) = 1 + 4tt 



n p f 
Jo 



h(r) 



sinirq) 2 



r dr 



rq 



(9) 



where n p is the number density, r is the center-to-center 
distance, q is the magnitude of the scattering vector given 
by q = 4>7i; X~ l sin(6 /2) (where 6 is the total scattering 
angle) and h(r) is the total correlation function which is 
given by h(r) = g(r) — 1. The radial distribution function 
was computed from BD simulations using the center-to- 
center protein distances. We estimated the convergence 
of the g(r) by checking that it was not varying with 
increasing simulation time. This was done by computing 
the g(r) over the full trajectory and comparing this g(r) 
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with an average g(r) computed from 20 segments selected 
sequentially from the trajectory. 

Test systems of two spherical particles 

For a system composed of two charged soft-sphere par- 
ticles interacting via a Debye-Huckel potential, the long- 
range contribution to the second virial coefficient can 
be computed by integrating Equation 6. This equation 
can be solved analytically by expanding the exponen- 
tial e~ w(x ^l kBT up to the second order and substituting 
the Debye-Huckel expression for the potential of mean 
force [29,34]. 

Only the long-range contribution to the second virial 
coefficient is taken into account in the analysis. Hence, the 
lower bound of the integration {lb) is not 0 but it is set to 
the sum of the protein radii (di+aj) plus one or two Debye 
lengths (1/k). For example, solving Equation 5 setting the 
lower bound to lb = {ai + aj) + 1/k gives 



1/k _ Z i Z ) e l 



B 22 - 



2ep 



2 + KCL K 
1 + KCL 47T€o€ r ekT (1 + KO) 2 



ZiZje t 



(10) 



where e is the base of the natural logarithm, e/ is the ele- 
mentary charge and p is the concentration of the ions 
(equivalent to the ionic strength for monovalent ions). 

The reason for considering only the long-range con- 
tribution is two-fold. Firstly, our purpose is to assess 
the accuracy of the long-range Debye-Huckel potential 
included in the BD simulation model. Secondly, for the 
expansion of the exponential e~ w l kT up to the second 
order to be reasonably accurate, \w/kT\ <$C 1 is required. 
This means that the short-range contribution of B22 at 
low ionic strength or for highly charged systems cannot be 
obtained using Equation 5. 

In the numerical integration, the two particles were 
represented by spherical fullerene-like particles of radius 
6 A composed of 180 atoms. A partial point charge was 
placed on each atom. The total charge of each sphere 
was uniformly distributed over all the atoms. Different 
systems were simulated by varying the net charge and 



the ionic strength (see Table 1 and Table 2 in Results 
and discussion). The interaction energy between the two 
particles is given by 



h 

+ \ ^2 ° e/ 2 ( r ;i) ' ®i [electrostatic] 



ZiZiefe-^-^ 

_| J —L 

47T€o€ r r(l + KCL) 

~l~ ^ ] Esoftcorei ( r W2 ) 



[Debye-Huckel ] 



(ID 



+ ^ E softcore 2 ( r «i ) [soft-core repulsion] 

To compute the second virial coefficient, one particle 
was kept fixed at the center of the simulation box and 
the other was moved on a regular lattice within the sim- 
ulation box, avoiding overlaps with the central particle. 
The size of the box was set to 400 x 400 x 400 A 3 and 
the dimension of the lattice was set to 100 x 100 x 100 
vertices. The interaction energy (Equation 11) was com- 
puted for each position assumed by the second particle 
and the second virial coefficient was computed by inte- 
grating Equation 6 numerically with the potential of mean 
force, w(r) = AG]^ , where r is the center- to-center 
distance. As for the analytical computation of B22> the 
integration was performed setting half, one, or two Debye 
lengths as the lower bound of the integral. 

We considered two spherical particles i and j with cor- 
responding radii at and aj and net charges Zi and Zj, each 
resulting from 180 partial point charges uniformly dis- 
tributed near the surface of each particle at a distance r 
from the particles center. Six different combinations of 
net charges on the particles were tested, namely: +1/+1, 
+5/+5, +10/+10 and +1/-1, +5/-5, +10/-10 (in units 
of elementary charge). For each pair of particles the inte- 
gration was performed at different ionic strengths, 5 mM 
and 300 mM. These two ionic strengths were chosen 



Table 1 Long range contribution to the £22 value at 5 mM ionic strength for the two soft-sphere systems 









) 


r (EP 2 ooa) 




3 22 




R (A) 
B 22 




(e/) 


(e/) 


1/JC 


2/* 


1/JC 


2/* 


V* 


2/K 


1/JC 


2/k 


+ 1 


+ 1 


0.0 


0.0 


] 3.0(0-0) 


0.0 


29.0(0-°) 


14.1(0-°) 


30.4 


16.1 


+ 1 


-1 


0.0 


0.0 


_ 133(0.0) 


0.0 


-29.4 (0 - 0) 


-14.2(0-°) 


-30.7 


-16.2 


+5 


+5 


0.0 


0.0 


261.5< ai > 


0.0 


636.9(0.1) 


339.3(0.2) 


675.8 


397.4 


+5 


-5 


0.0 


0.0 


-432.7((B) 


0.0 


-853.8(° 3 ) 


_367.8(°- 3 ) 


-878.0 


-429.8 


+10 


+ 10 


0.0 


0.0 


603.6(°- 2 ) 


0.0 


1916.9 (a2) 


1216.5< a7 > 


1490.0 


1428.1 


+10 


-10 


0.0 


0.0 


-5037.6< 9 - 7 ) 


0.0 


_7 338 .9(10.5) 


1682.2( 1 - 7 ) 


-4738.0 


-1867.6 



Computed values with electrostatic potential grids only (EP) are compared to the computed values with electrostatic potential grids plus the Debye-Huckel potential 
(EP+DHP) and to the analytical values (A). Lower bounds to integrate the B 2 2 were assigned as one or two Debye lengths see manuscript for details. The B 2 2 
values are given in units of (mol mL g~ 2 ) x 1 0 4 with standard deviations in parentheses. 



Mereghetti etal. BMC Biophysics 2014, 7:4 
http://www.biomedcentral.eom/2046-1682/7/4 



Page 6 of 1 1 



Table 2 Long range contribution to the B 2 2 values at 300 mM ionic strength for the two soft-sphere systems 
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B 22 
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B 22 




B 22 


ooa+DHP) 
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2/K 


I/* 


2/K 


I/* 


2/K 


V* 


2/K 


+ 1 


+ 1 


0.3< ao > 


0.1(0.0) 


03 (0.0) 


0.1(0.0) 


0.3(0-0) 


0.1(0.0) 


0.3 


0.1 


+ 1 


-1 


-0.3< ao > 


_ 0.1 (0-0) 


_0.3< ao > 


0.1(0.0) 


0.3(0-0) 


01 (o.o) 


-0.3 


-0.1 


+5 


+5 


6.4(0.0) 


33 (0.0) 


6 . 6 (0.D 


3.5(0.0) 


6.7(0-0) 


3.6(0-0) 


7.4 


4.0 


+5 


-5 


_ 9-4 (o.o) 


-3.7(0.0) 


_9. 6 (02) 


-3.9(0-1) 


-9.7(0.1) 


4.0(0.D 


-11.6 


-4.6 


+ 10 


+ 10 


181 (o.o) 


1 1 .4 (ao) 


1 8 . 8 (0.1) 


12. 0 (°- 1 ) 


19.4(0-1) 


12.6(°- 1 ) 


4.2 


12.8 


+ 10 


-10 


_ 97 . 6 (0.7) 


_ 18.4(0.1) 


_ 98 .1(0.4) 


_ 19.0(0-1) 


-98.9(° 7 ) 


_ 19.6(0.1) 


-72.4 


-22.0 



Computed values with electrostatic potential grids only (EP) are compared to the computed values with electrostatic potential grids plus the Debye-Huckel potential 
(EP+DHP) and to the analytical values (A). Lower bounds to integrate the B 2 i were set to one or two Debye lengths see manuscript for details. B 2 2 values are 
given in unit of (mol ml_ g~ 2 ) x 1 0 4 with standard deviations in parentheses. 



to assess the importance of the Debye-Huckel term at 
low and high salt conditions (compared to the 150 mM 
physiological ionic strength). The computed values were 
obtained by with and without inclusion of the Debye- 
Huckel potential 

From the set of approximately 10 6 interaction energies 
computed at the lattice vertices (avoiding the overlapping 
region), we extracted 100 random subsets of 10 5 values. 
For each subset, the second virial coefficient was com- 
puted. Then, an average .622 and a standard deviation over 
the subset was calculated. 

BD Simulations of protein solutions 

BD simulations were performed with SDAMM [3], a par- 
allelized program based on the SDA software [8] capable 
of handling many proteins (10 3 -10 4 ) treated as rigid bod- 
ies in atomic detail. For further details, see [3]. 

BD simulations were carried out for 250 protein 
molecules that were initially randomly positioned (avoid- 
ing overlaps) in a cubic box with periodic boundary con- 
ditions. The dimensions of the simulation box were varied 
according to the concentration of the protein solution. 

The Debye-Huckel interaction between a pair of pro- 
teins was computed up to a distance cutoff of 4 times 
the side of the electrostatic grid. If the simulation box 
was small, to avoid self-image interactions, this cutoff was 
rescaled to a value equal to half of the simulation box size. 

Each system was subjected to 5 or 10 fis of simulation 
at 300 K. Equilibration was assessed by monitoring the 
convergence of the radial distribution function and the 
stabilization of the energies. In all cases, 1 /is was suffi- 
cient to obtain an equilibrated system according to these 
criteria and the remaining 4 or 9/xs were used for the anal- 
ysis. The integration timestep was 0.5 ps. The positions 
and orientations of the proteins were recorded along with 
energy values every 0.5 ns. 

Simulations of HEWL were performed at 14, 28, 57 
and 85 g/L for comparison with experimental long-time 
translational self-diffusion coefficients [35]. Four sets of 



simulations were performed varying the ionic strength 
(1 mM and 5 mM) and including or omitting the analytical 
Debye-Huckel potential. Simulations were performed for 

5 /is. 

Simulations of BSA were performed at 0.9, 4.5, 9, 
18, 45, 90 g/L for comparison with the experimental 
small angle X-ray scattering (SAXS) intensities described 
in ref. [32]. Two sets of simulations were performed. 
In one set, the Debye-Huckel potential was included, 
whereas in the other set, the Debye-Huckel poten- 
tial was omitted. Because of the faster convergence of 
the higher concentration simulations, simulations at 0.9, 
4.5, 9 and 18 g/L were performed for 10 /is whereas 
the simulations at 45 and 90 g/L were performed for 
5 /is. 

Protein preparation 

The crystal structure of hen egg white lysozyme (HEWL) 
was taken from the Protein Data Bank (ref): lhel. The 
structure of BSA used for the simulations was a model 
taken from Modbase [36]. It was obtained by homology 
modelling based on the crystal structure of human serum 
albumin (HS A) [37]. 

Polar hydrogen atoms were added to the structures 
according to the specified pH and ionic strength (IS) using 
the H++ software [38]. The simulations of HEWL were 
performed at pH 5 ; the computed net charge of HEWL 
was +10e. The simulations of BSA were performed at pH 
7. BSA had a computed net charge of -16e. 

Atomic partial charges and radii were assigned to all 
the atoms from the OPLS united atom force field [39]. 
Electrostatic potential grids O were computed by solving 
the linearized Poisson-Boltzmann equation using the pro- 
gram UHBD [40]. The grid size was set to 100 x 100 x 
100 A 3 for HEWL and 200 x 200 x 200 A 3 for BSA with 
a grid spacing of 1.0 A. Non-polar desolvation, electro- 
static desolvation and soft-core repulsion grids were set to 
100 x 100 x 100 A 3 for HEWL and 130 x 130 x 130 A 3 for 
BSA, with a grid spacing of 1.0 A. 
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Results and discussion 

Comparison of simulations and analytical results for 
systems of two spherical particles 

The two spheres system (see Computational Details 
section) was simulated with different combinations of 
net solute charge at two ionic strengths with and with- 
out inclusion of the Debye-Huckel potential For each 
system, the analytical value of the long range con- 
tribution to the ,622 was compared to the computed 
one. All values are given in Table 1 for 5 mM and 
Table 2 for 300 mM ionic strength. For a better com- 
prehension of the length scale of the contribution of 
the electrostatic potential to the second virial coef- 
ficient, the analytical .622 values from the analytical 
calculations and from the simulations were obtained 
using different lower bounds for integrating Equation 6. 
We first consider the systems at low ionic strength 
(5 mM). 

5 mM ionic strength 

Let us first consider the integration done with a lower 
bound of one Debye length which at 5 mM ionic strength 
corresponds to 43 A. From Table 1, it is clear that when 
using a grid of 100 x 100 x 100 A 3 without the Debye- 
Huckel potential, the long-range decay of the electrostatic 
potential is not captured. This result is expected since the 
electrostatic potential grid size is of the same order as the 
Debye length. Doubling the length of the side of the grid 
results in a B22 value that is approximately 50% of the ana- 
lytical value. The long-range tail (beyond 100 A) of the 
electrostatic potential is missing and it is apparent that it 
represents an important contribution to the second virial 
coefficient. 

By turning on the Debye-Huckel potential and keeping 
the smaller electrostatic potential grid (side length: 100 A), 
more than the 90% of the analytical B22 value is recov- 
ered. For systems with the highest net charge at one Debye 
length, the potential is too high and the integral expression 
in Equation 6 diverges. 

For a perfectly isotropic case, such as this one, the 
Debye-Huckel potential smoothly recovers the truncation 
of the electrostatic potential due to the finite grid. This can 
be seen from the electrostatic potential energy computed 
by varying the inter-particle separation (see Additional 
filel). 

At two Debye lengths (2//c), the ,622 value of the sys- 
tems with the smaller grid (100 A) without the Debye- 
Huckel potential is zero, since the grid is smaller than 
the Debye length. By doubling the grid dimension, the 
side of the grid becomes of the same order as the Debye 
length and the B22 is still not computed correctly. With 
the Debye-Huckel potential and the smaller grid, however, 
the analytical second virial coefficient can be reproduced 
well. 



300 mM ionic strength 

Increasing the ionic strength to 300 mM, at lower bounds 
of one or two Debye lengths (5.5 A), the B22 values com- 
puted using only the smaller electrostatic potential grid 
agree rather well with the analytical values, see Table 2. 
Doubling the grid dimensions or adding the Debye- 
Huckel potential is not required because more than 90% 
of the interactions are captured within one Debye length. 
It is clear that at 300 mM ionic strength, the grid-based 
formalism is sufficient to properly describe the long-range 
electrostatic interaction, even using the smaller grid. 

Protein systems modeled in atomic detail 

We now turn to more complex and realistic systems 
composed of solutions of proteins represented in atomic 
detail subjected to BD simulation as described in the 
Computational Details section. 

Scattering intensities Several BSA solutions at different 
concentrations were simulated for 10 /xs to 20 /xs using 
BD. To assess the effect of the Debye-Huckel approxima- 
tion on the BSA self-interactions, two sets of simulations 
were performed. In one set, the Debye-Huckel potential 
was included whereas in the other set, it was omitted. 

Normalized small angle scattering intensities were com- 
puted using Equation 8 and compared to experimental 
SAXS intensities. The experiments were performed with- 
out added salt which corresponds to an ionic strength up 
to 5 mM [31,32]. This non-zero ionic strength arises from 
several factors such as dissolved CO2, a residual amount 
of salt present in the protein solution, and the dissociation 
of surface groups upon solvation [31,32]. Simulations were 
performed at 5 mM ionic strength with a corresponding 
Debye length of 43.1 A. 

As shown in Figure 1, the scattering intensities obtained 
from the simulations with the Debye-Huckel approxima- 
tion reproduce experimental SAXS intensities better then 
the intensities calculated from simulations which do not 
include the Debye-Huckel interaction. In particular, the 
greatest improvement is seen at low q values, i.e. long 
range interactions are accurately captured. At high con- 
centrations, the Debye-Huckel approximation tends to 
overestimate the height of the correlation peak seen in the 
normalized experimental intensities. This phenomenon 
can be explained considering that simulations have been 
performed at 5 mM ionic strength, but that at high protein 
concentrations, the effective ionic strength may be higher 
due to the presence of highly charged proteins. Indeed, 
the correlation peak is lower in the simulations without 
the Debye-Huckel approximation (see also Figure 2 and 
Figure 3). This suggests that at low ionic strength and 
high protein concentration, the ionic strength of the sim- 
ulation should be slightly increased to better reproduce 
experimentally observed scattering intensities. 
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Figure 1 BSA SAS intensities. Experimental [32] (dashed lines) and computed (continuous lines) normalized small angle scattering intensities at 
different concentrations (indicated on the plots) of BSA. Computed curves from simulations without (A) and with (B) the Debye-Huckel 
approximation. Curves are shifted by 0.2 on the vertical axis for better visibility. 



The computed static structure factors obtained from 
the two sets of simulations are compared in Figure 2. 
Focusing on the low q region (q < 0.1 nm~ l ), for a given 
concentration, the value of S(q) is lower when the Debye- 
Huckel potential is used. The long wavelength limit of 
S(q) is proportional to the normalized isothermal osmotic 
compressibility, vis.: 

lim S(q) = n p k B TxT 

where xt is the isothermal osmotic compressibility. 
(In the canonical ensemble, xt = — ^(fn)r = 

\ n P (f>?~) } ^ n P * s ^ e P rotem numDer density, and 
I<b is the Boltzmann constant [32,41,42]. The decrease of 



S(q) at low q values can be explained by the decrease 
of the osmotic compressibility due to the long-range 
electrostatic repulsion introduced with the Debye-Huckel 
potential [43]. 

The first peak in the S(q) represents the correlation 
between a pair of proteins. We observe that the simula- 
tions which include the Debye-Huckel potential show a 
shift of the first peak to lower q values (at high concen- 
trations) or the appearance of a peak (at low concentra- 
tions), indicating the presence of a long-range correlation 
between the proteins. With increasing concentration, the 
peak shifts to higher q values, suggesting a decrease of the 
correlation distance. The same effect can be seen better 
in real space from the radial distribution functions plotted 




0.2 0.4 0.6 0.8 1 
q (nm" ) 

Figure 2 BSA structure factors. Experimental [32] (dashed lines) and 
computed (continuous lines) structure factors at various 
concentrations (indicated on the plot) of BSA obtained from 
simulations without (dark green) and with (dark red) the 
Debye-Huckel approximation. Curves are shifted by 0.2 on the vertical 
axis for better visibility. 




Figure 3 BSA radial distribution functions. Computed radial 
distribution functions at various concentrations (indicated on the plot) 
of BSA obtained from simulations without (dark green) and with (dark 
red) the Debye-Huckel approximation. Curves are shifted by 0.2 on 
the vertical axis for better visibility. Averages and standard deviations 
of the g(r) are shown by the dark line and light color, respectively. 
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in Figure 3 where it can be seen that the introduction of a 
long-range repulsion pushes the proteins away from each 
other. It also leads to a more structured solution, with the 
appearance of a second peak in the simulations at 90 g/L 
protein concentration. 

Long-time self-diffusion coefficients Besides the effect 
on protein-protein interactions, the addition of the 
Debye-Huckel potential also has consequences for the 
dynamics of the proteins. Simulations of HEWL were per- 
formed at low ionic strength (1 and 5 mM) at different 
lysozyme concentrations and compared to experimen- 
tal diffusion coefficients obtained from pulsed gradient 
spin echo NMR for HEWL solutions without added 
salt at pH 4.9. As shown in Figure 4, the presence 
of the Debye-Huckel potential systematically lowers the 
long-time self-diffusion coefficients. This effect can be 
explained considering that, for a given concentration, sim- 
ulations which include the Debye-Huckel potential cor- 
respond to a larger effective concentration due to the 
long-range repulsive interaction [43,44], In general, the 
magnitude of the effect on the diffusion coefficient due to 
the Debye-Huckel potential is related to the ionic strength 
of the solution, the size of the protein, and the protein 
concentration. For proteins whose size is comparable to 
the Debye length, as in our case, this effect can be 
significant. For very large proteins, the Debye length can 
be much smaller than the size of the protein, and hence, 
adding the long-range Debye-Huckel interaction may lead 
only to small effects on the diffusion coefficient. 

Simulations performed at 1 mM ionic strength under- 
estimate the diffusion coefficients compared to the exper- 
imental values (see Figure 4). As described above for the 
BSA case, the ionic strength of the solution is affected 
by several factors. Thus, it is possible that the value of 



1 mM used in the simulations does not correctly describe 
the effective ionic strength of the experimental solutions. 
We therefore also performed simulations at higher ionic 
strength (5 mM), obtaining better agreement with the 
experimental data, see Figure 4. 

Methodological considerations 

The Debye-Huckel potential has been implemented 
together with cubic grids for the proteins. The transition 
from the gridded potential to the Debye-Huckel poten- 
tial with increasing distance from a solute center occurs 
at the shortest distance to the grid boundary. Thus, cubic 
grids permit the most efficient implementation of the 
Debye-Huckel correction. Their use is usually appropri- 
ate for globular proteins, however, it may become an issue 
when modeling large elongated molecules. For the lat- 
ter, a large number of grid points on a cubic grid will 
have very low (negligible) values of the mapped interac- 
tion potentials, leading to an unnecessarily high memory 
requirement. 

On the other hand, an advantage of the Debye-Huckel 
implementation is that it removes the requirement for the 
electrostatic potential to have very small values at the grid 
edges; the electrostatic potential is only required to be 
centrosymmetric. This means that smaller grids can be 
used with the long-range interactions being captured by 
the Debye-Huckel with only a small computational cost 
(see Additional file 2). 

Using the Debye-Huckel correction may be an issue for 
some highly or non-uniformly charged systems as it can 
lead to force discontinuities at the grid boundaries. A pos- 
sible solution to this problem, currently not implemented, 
is to apply an interpolation function between the electro- 
static potential grid and the Debye-Huckel potential for 
computing the forces at the grid boundary. 
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Figure 4 HEWL translational diffusion coefficients. Normalized long-time translational self-diffusion coefficients of HEWL at low ionic strength. 
Simulations were performed at 1 mM (A) and 5 mM (B) ionic strength. Experimental values from ref. [35] (black diamonds), and computed values 
from BD simulations with (red squares) and without (green squares) Debye-Huckel potential are shown. The Tokuyama [22] analytical model is 
shown by the black dotted line. Insets are log-log plots of the same data. 
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Conclusions 

We have here described the implementation of a Debye- 
Huckel correction for the computation of grid-based 
electrostatic interaction energies and forces for use in 
atomically detailed many-protein Brownian dynamics 
simulations. The ability of this many-protein BD method 
to correctly reproduce small angle scattering data and 
diffusion coefficients, was previously shown for several 
proteins [3,12]. Due to computational limitations on the 
size of the electrostatic interaction grids, the method 
could not be applied to highly charged systems or low 
ionic strength conditions without impairing the accuracy 
of the resulting simulations. The introduction of the sim- 
ple Debye-Huckel correction described in this paper with 
its very low associated computational costs allowed us to 
extend the range applicability of this BD method to highly 
charged systems at low ionic strength. In particular, com- 
parison of the model with the Debye-Huckel correction to 
analytical results for spherical solutes, as well as to exper- 
imental SAXS intensities for BSA protein solutions, and 
to long-time self-diffusion coefficients of HEWL protein 
solutions, showed good agreement. Some other potential 
applications of the methodology are the simulation of pro- 
tein crystallization, of protein-surface adsorption, and of 
heterogeneous crowded protein solutions. Furthermore, 
the Debye-Huckel correction described here should be of 
value in implicit solvent molecular dynamics simulations 
which make use of gridded interaction potentials [13-16]. 

Additional files 



Additional file 1 : Charged spheres electrostatic potential energy. 

Electrostatic potential energy of two uniformly charged spheres for 
different net charge (e) combinations. Panels A,B: +1/-1; panels C,D: 

+5/+5, +5/-5; panels E,F: +1 0/+1 0, +1 0/-1 0. Red and green colors show 
interactions between charges of the same and opposite sign, respectively. 
The Debye-Huckel analytical approximation, Equation 4 (continuous line), 
Brownian dynamics without Debye-Huckel term (crosses) and Brownian 
dynamics with Debye-Huckel term (circles) are shown. Interactions are 
computed at 5 mM (left panels) and 300 mM (right panels). The abscissa is 
the particle center-to-center distance divided by the particle diameter. The 
inclusion of the Debye-Huckel potential recovers the dependence of the 
energy on particle separation computed with the analytical model. 

Additional file 2: Runtime versus protein concentration for the 
simulations of BSA. The Debye-Huckel correction requires very little 
additional computational effort. Indeed, at low concentrations, the 
inclusion of the Debye-Huckel correction keeps the like-charged molecules 
apart and therefore reduces the number of pairs of molecules for which 
grid-type potentials are used to compute intermolecular forces, thus 
leading to reduced run-times. 
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